/**********************************************************************

  Audacity: A Digital Audio Editor

  FreqWindow.cpp

  Dominic Mazzoni

  *******************************************************************//**

  \class FreqWindow
  \brief Displays a spectrum plot of the waveform.  Has options for
  selecting parameters of the plot.

  Has a feature that finds peaks and reports their value as you move
  the mouse around.

  *//****************************************************************//**

  \class FreqPlot
  \brief Works with FreqWindow to dsplay a spectrum plot of the waveform.
  This class actually does the graph display.

  Has a feature that finds peaks and reports their value as you move
  the mouse around.

  *//*******************************************************************/

/*
  Salvo Ventura - November 2006
  Extended range check for additional FFT windows
  */

#ifdef __BORLANDC__
#pragma hdrstop
#endif

#include "stdafx.h"
#include <math.h>
#include <string.h>
#include <memory.h>
#include "FFT.h"
#include "FreqAlg.h"

#ifndef min
#define min(a,b)            (((a) < (b)) ? (a) : (b))
#endif

// FreqWindow

FreqAlg::FreqAlg()
: mData(NULL)
, mProcessed(NULL)
, mWindowSize(1024)
{
	mRate = 0;
	mDataLen = 0;

	mAlg = 0;
	mFunc = 0;

	dBRange = 120;
	mLogAxis = false;
}

FreqAlg::~FreqAlg()
{
	if (mData)
		delete[] mData;

	if (mProcessed)
		delete[] mProcessed;

	DeinitFFT();
}

void FreqAlg::SetRate(double rate)
{
	mRate = rate;
}

void FreqAlg::SetAlg(int alg)
{
	mAlg = alg;
}

void FreqAlg::SetFunc(int func)
{
	mFunc = func;
}

void FreqAlg::SetWindowSize(int size)
{
	mWindowSize = size;
}

void FreqAlg::SetAxis(int axis)
{
	mLogAxis = axis ? true : false;
}

void FreqAlg::SetData(float* data, int len)
{
	if (mData) {
		delete[] mData;
	}

	mDataLen = len;
	mData = new float[mDataLen];

	memset(mData, 0, mDataLen * sizeof(float));
	memcpy(mData, data, mDataLen * sizeof(float));
}

void FreqAlg::SetData(short* data, int len)
{
	double fscale = 1.0 / 32768.0;
	float* fdata = new float[len];
	for (int i = 0; i < len; i++)
	{
		fdata[i] = (float)((double)data[i] * fscale);
	}

	SetData(fdata, len);
	delete[] fdata;
}

float FreqAlg::CubicInterpolate(float y0, float y1, float y2, float y3, float x)
{
	float a, b, c, d;

	a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
	b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
	c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
	d = y0;

	float xx = x * x;
	float xxx = xx * x;

	return (a * xxx + b * xx + c * x + d);
}

float FreqAlg::CubicMaximize(float y0, float y1, float y2, float y3, float * max)
{
	// Find coefficients of cubic

	float a, b, c, d;

	a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
	b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
	c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
	d = y0;

	// Take derivative

	float da, db, dc;

	da = 3 * a;
	db = 2 * b;
	dc = c;

	// Find zeroes of derivative using quadratic equation

	float discriminant = db * db - 4 * da * dc;
	if (discriminant < 0.0)
		return float(-1.0);              // error

	float x1 = (-db + sqrt(discriminant)) / (2 * da);
	float x2 = (-db - sqrt(discriminant)) / (2 * da);

	// The one which corresponds to a local _maximum_ in the
	// cubic is the one we want - the one with a negative
	// second derivative

	float dda = 2 * da;
	float ddb = db;

	if (dda * x1 + ddb < 0)
	{
		*max = a*x1*x1*x1 + b*x1*x1 + c*x1 + d;
		return x1;
	}
	else
	{
		*max = a*x2*x2*x2 + b*x2*x2 + c*x2 + d;
		return x2;
	}
}

float FreqAlg::GetProcessedValue(float freq0, float freq1)
{
	int alg = mAlg;

	float bin0, bin1, binwidth;

	if (alg == 0) {
		bin0 = freq0 * mWindowSize / mRate;
		bin1 = freq1 * mWindowSize / mRate;
	}
	else {
		bin0 = freq0 * mRate;
		bin1 = freq1 * mRate;
	}
	binwidth = bin1 - bin0;

	float value = float(0.0);

	if (binwidth < 1.0) {
		float binmid = (bin0 + bin1) / 2.0;
		int ibin = int(binmid) - 1;
		if (ibin < 1)
			ibin = 1;
		if (ibin >= mProcessedSize - 3)
			ibin = mProcessedSize - 4;

		value = CubicInterpolate(mProcessed[ibin],
			mProcessed[ibin + 1],
			mProcessed[ibin + 2],
			mProcessed[ibin + 3], binmid - ibin);

	}
	else {
		if (int(bin1) > int(bin0))
			value += mProcessed[int(bin0)] * (int(bin0) + 1 - bin0);
		bin0 = 1 + int(bin0);
		while (bin0 < int(bin1)) {
			value += mProcessed[int(bin0)];
			bin0 += 1.0;
		}
		value += mProcessed[int(bin1)] * (bin1 - int(bin1));

		value /= binwidth;
	}

	return value;
}

void FreqAlg::Recalc()
{
	if (mProcessed)
		delete[] mProcessed;
	mProcessed = NULL;

	if (!mData) {
		return;
	}

	int f = 10;
	int alg = mAlg;
	int windowFunc = mFunc;
	long windowSize = mWindowSize;

	if (!(windowSize >= 32 && windowSize <= 65536 &&
		alg >= 0 && alg <= 4 && windowFunc >= 0 && windowFunc < f)) {
		return;
	}

	if (mDataLen < mWindowSize) {
		return;
	}

	mProcessed = new float[mWindowSize];

	int i;
	for (i = 0; i < mWindowSize; i++)
		mProcessed[i] = float(0.0);
	int half = mWindowSize / 2;

	float *in = new float[mWindowSize];
	float *in2 = new float[mWindowSize];
	float *out = new float[mWindowSize];
	float *out2 = new float[mWindowSize];
	float *win = new float[mWindowSize];

	// initialize the window
	for (int i = 0; i < mWindowSize; i++)
		win[i] = 1.0;
	WindowFunc(windowFunc, mWindowSize, win);
	// Scale window such that an amplitude of 1.0 in the time domain
	// shows an amplitude of 0dB in the frequency domain
	double wss = 0;
	for (int i = 0; i<mWindowSize; i++)
		wss += win[i];
	if (wss > 0)
		wss = 4.0 / (wss*wss);
	else
		wss = 1.0;

	int start = 0;
	int windows = 0;
	while (start + mWindowSize <= mDataLen) {
		for (i = 0; i < mWindowSize; i++)
			in[i] = win[i] * mData[start + i];

		switch (alg) {
		case 0:                  // Spectrum
			PowerSpectrum(mWindowSize, in, out);

			for (i = 0; i < half; i++)
				mProcessed[i] += out[i];
			break;

		case 1:
		case 2:
		case 3:   // Autocorrelation, Cuberoot AC or Enhanced AC

			// Take FFT
#ifdef EXPERIMENTAL_USE_REALFFTF
			RealFFT(mWindowSize, in, out, out2);
#else
			FFT(mWindowSize, false, in, NULL, out, out2);
#endif
			// Compute power
			for (i = 0; i < mWindowSize; i++)
				in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);

			if (alg == 1) {
				for (i = 0; i < mWindowSize; i++)
					in[i] = sqrt(in[i]);
			}
			if (alg == 2 || alg == 3) {
				// Tolonen and Karjalainen recommend taking the cube root
				// of the power, instead of the square root

				for (i = 0; i < mWindowSize; i++)
					in[i] = pow(in[i], 1.0f / 3.0f);
			}
			// Take FFT
#ifdef EXPERIMENTAL_USE_REALFFTF
			RealFFT(mWindowSize, in, out, out2);
#else
			FFT(mWindowSize, false, in, NULL, out, out2);
#endif

			// Take real part of result
			for (i = 0; i < half; i++)
				mProcessed[i] += out[i];
			break;

		case 4:                  // Cepstrum
#ifdef EXPERIMENTAL_USE_REALFFTF
			RealFFT(mWindowSize, in, out, out2);
#else
			FFT(mWindowSize, false, in, NULL, out, out2);
#endif

			// Compute log power
			// Set a sane lower limit assuming maximum time amplitude of 1.0
			float power;
			float minpower = 1e-20*mWindowSize*mWindowSize;
			for (i = 0; i < mWindowSize; i++)
			{
				power = (out[i] * out[i]) + (out2[i] * out2[i]);
				if (power < minpower)
					in[i] = log(minpower);
				else
					in[i] = log(power);
			}
			// Take IFFT
#ifdef EXPERIMENTAL_USE_REALFFTF
			InverseRealFFT(mWindowSize, in, NULL, out);
#else
			FFT(mWindowSize, true, in, NULL, out, out2);
#endif

			// Take real part of result
			for (i = 0; i < half; i++)
				mProcessed[i] += out[i];

			break;
		}                         //switch

		start += half;
		windows++;
	}

	switch (alg) {
		double scale;
	case 0:                     // Spectrum
		// Convert to decibels
		mYMin = 1000000.;
		mYMax = -1000000.;
		scale = wss / (double)windows;
		for (i = 0; i < half; i++)
		{
			mProcessed[i] = 10 * log10(mProcessed[i] * scale);
			if (mProcessed[i] > mYMax)
				mYMax = mProcessed[i];
			else if (mProcessed[i] < mYMin)
				mYMin = mProcessed[i];
		}
		if (mYMin < -dBRange)
			mYMin = -dBRange;
		if (mYMax <= -dBRange)
			mYMax = -dBRange + 10.; // it's all out of range, but show a scale.
		else
			mYMax += .5;

		mProcessedSize = half;
		mYStep = 10;
		break;

	case 1:                     // Standard Autocorrelation
	case 2:                     // Cuberoot Autocorrelation
		for (i = 0; i < half; i++)
			mProcessed[i] = mProcessed[i] / windows;

		// Find min/max
		mYMin = mProcessed[0];
		mYMax = mProcessed[0];
		for (i = 1; i < half; i++)
		if (mProcessed[i] > mYMax)
			mYMax = mProcessed[i];
		else if (mProcessed[i] < mYMin)
			mYMin = mProcessed[i];

		mYStep = 1;

		mProcessedSize = half;
		break;

	case 3:                     // Enhanced Autocorrelation
		for (i = 0; i < half; i++)
			mProcessed[i] = mProcessed[i] / windows;

		// Peak Pruning as described by Tolonen and Karjalainen, 2000

		// Clip at zero, copy to temp array
		for (i = 0; i < half; i++) {
			if (mProcessed[i] < 0.0)
				mProcessed[i] = float(0.0);
			out[i] = mProcessed[i];
		}

		// Subtract a time-doubled signal (linearly interp.) from the original
		// (clipped) signal
		for (i = 0; i < half; i++)
		if ((i % 2) == 0)
			mProcessed[i] -= out[i / 2];
		else
			mProcessed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);

		// Clip at zero again
		for (i = 0; i < half; i++)
		if (mProcessed[i] < 0.0)
			mProcessed[i] = float(0.0);

		// Find new min/max
		mYMin = mProcessed[0];
		mYMax = mProcessed[0];
		for (i = 1; i < half; i++)
		if (mProcessed[i] > mYMax)
			mYMax = mProcessed[i];
		else if (mProcessed[i] < mYMin)
			mYMin = mProcessed[i];

		mYStep = 1;

		mProcessedSize = half;
		break;

	case 4:                     // Cepstrum
		for (i = 0; i < half; i++)
			mProcessed[i] = mProcessed[i] / windows;

		// Find min/max, ignoring first and last few values
		int ignore = 4;
		mYMin = mProcessed[ignore];
		mYMax = mProcessed[ignore];
		for (i = ignore + 1; i < half - ignore; i++)
		if (mProcessed[i] > mYMax)
			mYMax = mProcessed[i];
		else if (mProcessed[i] < mYMin)
			mYMin = mProcessed[i];

		mYStep = 1;

		mProcessedSize = half;
		break;
	}

	delete[]in;
	delete[]in2;
	delete[]out;
	delete[]out2;
	delete[]win;
}

int FreqAlg::GetFreqDbSize() const
{
	return mProcessedSize;
}

const float* FreqAlg::GetFreqDb() const
{
	return mProcessed;
}

void FreqAlg::GetFreqDb(float* buffer, int len)
{
	int nCount = min(mProcessedSize, len);
	if (nCount > 0)
	{
		memcpy(buffer, mProcessed, nCount * sizeof(float));
	}
}

float FreqAlg::GetYMin() const
{
	return mYMin;
}

float FreqAlg::GetYMax() const
{
	return mYMax;
}